rm(list=ls())
setwd('./Fig.1B_E')


library(vegan)
library(ggplot2)
require(stringi)
require(magrittr)
require(tibble)
require(agricolae)
library(ggpubr)
library(ggsignif)
require(cowplot)
require(patchwork)

otu_raw <- read.csv(file="lipid_mumian.csv" ,row.names=1)
head(otu_raw)
otu_raw%>%
  as.data.frame%>%
  set_rownames(NULL)%>%
  column_to_rownames(var="name")->otu_raw
#otu_raw <- read.table(file="otu.txt",sep="\t",header=T,check.names=FALSE ,row.names=1)
otu <- t(otu_raw)
otu[is.na(otu)] <- 0 
View(otu)
sapply(otu,class)

#otu<- round(otu,digits=0)
#View(otu1)
otu.distance <- vegdist(otu, method = 'bray')
df_nmds <- metaMDS(otu.distance, k = 2)
summary(df_nmds)

df_nmds_stress <- df_nmds$stress
df_nmds_stress
stressplot(df_nmds)

df_points <- as.data.frame(df_nmds$points)
df_points 
df_points$samples <- row.names(df_points)
df_points$samples
names(df_points)[1:2] <- c('NMDS1', 'NMDS2')
head(df_points)

p <- ggplot(df_points,aes(x=NMDS1, y=NMDS2))+
  geom_point(size=3)+
  theme_bw()
p


group <- read.table("group.txt", sep='\t', header=T)
group
colnames(group) <- c("samples","group")
#group$samples<-stri_replace_all(group$samples,"_",regex="-")
#group$group<-stri_extract(group$group,regex="^[A-Za-z]+_[0-9][a-z]")
dplyr::left_join(df_points,group,by="samples")->df
View(df)
head(df)



color=c("#983238","#d96458","#c56140","#f49a5a","#ebcf4e")
p1<-ggplot(data=df,aes(x=NMDS1,y=NMDS2))+
  theme_bw()+
  geom_point(aes(color = group), shape = 19, size=3)+
  theme(panel.grid = element_blank())+
  geom_vline(xintercept = 0,lty="dashed", size = 1, color = 'grey50')+
  geom_hline(yintercept = 0,lty="dashed", size = 1, color = 'grey50')+
  geom_text(aes(label=samples, y=NMDS2+0.03,x=NMDS1+0.03,
                vjust=0, color = group),size=3.5, show.legend = F)+
  stat_ellipse(data=df,
               geom = "polygon",level=0.9,
               linetype = 1,size=0.5,
               aes(fill=group),
               alpha=0.2)+
  scale_color_manual(values = color) +
  scale_fill_manual(values = color)+
  theme(axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12,angle=90),
        axis.text.y=element_text(size=10),
        axis.text.x=element_text(size=10),
        panel.grid=element_blank(),
        legend.position="none")+
  ggtitle(paste('Stress=',round(df_nmds_stress, 3)))

p1

p2 <- ggplot(df,aes(x=group,y=NMDS2))+
  stat_boxplot(geom = "errorbar", width=0.1,size=0.5)+
  geom_boxplot(aes(fill=group), 
               outlier.colour="white",size=0.5)+
  theme(panel.background =element_blank(), 
        axis.line=element_line(color = "white"),
        axis.text.y = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        legend.position = 'none')+
  xlab("") + ylab("")+
  scale_fill_manual(values=c("#983238","#d96458","#c56140","#f49a5a","#ebcf4e"))+
  geom_signif(comparisons = list(c("A","B"),
                                 c("A","C"),
                                 c("B","C")),
              map_signif_level = T, 
              test = t.test, 
              y_position = c(0.3,0.4,0.35),
              tip_length = c(c(0,0),
                             c(0,0),
                             c(0,0)),
              size=0.8,color="black")+
  theme(axis.title=element_blank())
p2

p3 <- ggplot(df,aes(x=group,y=NMDS1))+
  stat_boxplot(geom = "errorbar", width=0.1,size=0.5)+
  coord_flip()+
  geom_boxplot(aes(fill=group), 
               outlier.colour="white",size=0.5)+
  theme(panel.background =element_blank(), 
        axis.line=element_line(color = "white"),
        axis.text.x = element_blank(),
        panel.grid = element_blank(),
        axis.ticks = element_blank(),
        legend.position = 'none')+
  xlab("") + ylab("")+
  scale_fill_manual(values=c("#983238","#d96458","#c56140","#f49a5a","#ebcf4e"))+
  geom_signif(comparisons = list(c("A","B"),
                                 c("A","C"),
                                 c("B","C")),
              map_signif_level = T, 
              test = t.test, 
              y_position = c(0.48,0.62,0.55),
              tip_length = c(c(0,0),
                             c(0,0),
                             c(0,0)),
              size=0.8,color="black")
p3
#ggarrange(p3, NULL, p1+theme(legend.position="none"), p2, widths = c(5,2), heights = c(2,4), align = "hv")

####Shannon diversity####
diversity(otu,index="shannon")%>%
  as.data.frame%>%
  set_names("H")%>%
  rownames_to_column(var="group")%>%
  dplyr::mutate(treatment=group%>%
                  stri_extract(regex="^[A-Za-z]+"),
                repetition=group%>%
                  stri_split(regex="_")%>%
                  sapply(function(x)x[length(x)]))->shannonDiversity
model<-aov(H~treatment,data=shannonDiversity)
out<-HSD.test(model,"treatment",group=TRUE,console=TRUE)
out$groups%>%
  as.data.frame%>%
  rownames_to_column(var="treatment")%>%
  dplyr::left_join(shannonDiversity%>%
                     dplyr::group_by(treatment)%>%
                     dplyr::summarise(Hmax=max(H)),
                   by="treatment")->treatmentGroups
c(min(shannonDiversity$H),max(shannonDiversity$H)+0.1*max(shannonDiversity$H))->yLimits
ggplot(data=shannonDiversity,
       aes(x=treatment,
           y=H,
           colour=treatment))+
  geom_point(show.legend=FALSE)+
  geom_boxplot(show.legend=FALSE,
               alpha=0.3)+
  geom_text(data=treatmentGroups,
            aes(x=treatment,
                y=Hmax,
                label=groups),
            vjust=-1,
            show.legend=FALSE)+
  scale_colour_manual(values=c("#983238","#d96458","#c56140","#f49a5a","#ebcf4e"))+
  scale_y_continuous(limits=yLimits)+
  labs(x="",
       y="Shannon Index")+
  theme_classic()+
  theme(axis.title.x=element_blank())->shannonPlot

#plot_grid(p3+
#            theme(plot.margin=margin(0,-1,-1,-1,unit="cm")),
#          NULL,
#          p1+
#            theme(plot.margin=margin(-1,-1,0.1,0.1,unit="cm")),
#          p2+
#            theme(plot.margin=margin(-1,-1,-1,-1,unit="cm")),
#          shannonPlot,
#          rel_widths=c(5,2),
#          rel_heights=c(2,4,4),
#          ncol=2,
#          align="hv")

design<-"AAA##
         AAA##
         BBBCC
         BBBCC
         BBBCC
         DDDDD
         DDDDD
         DDDDD"

p3+
  theme(plot.margin=margin(0,-1,-1,-1,unit="cm"))+
  p1+
  theme(plot.margin=margin(-1,-1,0.1,0.1,unit="cm"))+
  p2+
  theme(plot.margin=margin(-1,-1,-1,-1,unit="cm"))+
  shannonPlot+
  theme(plot.margin=margin(-1,-1,0.1,0,unit="cm"))+
  plot_layout(design=design)+

  plot_annotation(tag_levels="A",
                  tag_prefix="(",
                  tag_suffix=")")&
  theme(plot.tag=element_text(size=unit(14,"pt"),family="Arial"))->combnPlots
combnPlots
ggsave("combn_plots_20240524.svg",combnPlots)
ggsave("combn_plots_20240524.png",combnPlots,width=7,height=7)




  